#clean up
rm(list=ls())

#preferelaes
options(scipen=8)

#load required libraries
library(spdep) 
library(maps)
library(foreign)
library(maptools)
library(spatstat)
library(mgcv)
library(ggmap)
library(geoR)
library(akima)
library(scatterplot3d)
library(sandwich) 
library(multiwayvcov) 
library(lmtest)
library(stringr)
library(xtable)
library(lattice)

####LOAD DATA###
#setwd("/Volumes/MONOGAN/registrationSpatial/LA/") ###SET TO YOUR WORKING DIRECTORY
la<-read.dta('laRegistration.dta')
la<-la[order(la$year,la$parish),]

###DESCRIPTIVES### 
summary(la)
boxplot(la$pctrep~la$year,xlab="Year",ylab="Proportion GOP Registration",ylim=c(0,.5))
cor(la$lnOil1,la$lnGas1)

###Choropleth map of GOP registration###
#front matter common to all maps
n.col<-5
br <- c(0,.05,.1,.2,.3,1)
shading<-gray((n.col-1):0/(n.col-1))

#county list common to all maps
counties.to.map<-str_trim(paste('louisiana', str_replace(tolower(la$parish[la$year==1966]),pattern="[.]",replacement=""), sep=','))
counties.to.map[50]<-"louisiana,st martin:north" 
counties.to.map[65]<-"louisiana,st martin:south" 
counties.to.map<-sort(counties.to.map)
#counties.to.map

#1965 map
X.1965 <- la$pctrep_t1[la$year==1966][c(1:49,50,50,51:64)]
X.1965.grp<-findInterval(X.1965, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1965.shad<-shading[X.1965.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1965.shad, exact=T)

#1966 map
X.1966 <- la$pctrep[la$year==1966][c(1:49,50,50,51:64)]
X.1966.grp<-findInterval(X.1966, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1966.shad<-shading[X.1966.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1966.shad, exact=T)

#1970 map
X.1970 <- la$pctrep[la$year==1970][c(1:49,50,50,51:64)]
X.1970.grp<-findInterval(X.1970, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1970.shad<-shading[X.1970.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1970.shad, exact=T)

#1976 map
X.1976 <- la$pctrep[la$year==1976][c(1:49,50,50,51:64)]
X.1976.grp<-findInterval(X.1976, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1976.shad<-shading[X.1976.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1976.shad, exact=T)

#1977 map
X.1977 <- la$pctrep[la$year==1977][c(1:49,50,50,51:64)]
X.1977.grp<-findInterval(X.1977, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1977.shad<-shading[X.1977.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1977.shad, exact=T)

#1978 map
X.1978 <- la$pctrep[la$year==1978][c(1:49,50,50,51:64)]
X.1978.grp<-findInterval(X.1978, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1978.shad<-shading[X.1978.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1978.shad, exact=T)

#1980 map
X.1980 <- la$pctrep[la$year==1980][c(1:49,50,50,51:64)]
X.1980.grp<-findInterval(X.1980, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1980.shad<-shading[X.1980.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1980.shad, exact=T)

#1986 map
X.1986 <- la$pctrep[la$year==1986][c(1:49,50,50,51:64)]
X.1986.grp<-findInterval(X.1986, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1986.shad<-shading[X.1986.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1986.shad, exact=T)

#1990 map
X.1990 <- la$pctrep[la$year==1990][c(1:49,50,50,51:64)]
X.1990.grp<-findInterval(X.1990, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1990.shad<-shading[X.1990.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1990.shad, exact=T)

#1996 map
X.1996 <- la$pctrep[la$year==1996][c(1:49,50,50,51:64)]
X.1996.grp<-findInterval(X.1996, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.1996.shad<-shading[X.1996.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1996.shad, exact=T)

#2000 map
X.2000 <- la$pctrep[la$year==2000][c(1:49,50,50,51:64)]
X.2000.grp<-findInterval(X.2000, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.2000.shad<-shading[X.2000.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2000.shad, exact=T)

#2006 map
X.2006 <- la$pctrep[la$year==2006][c(1:49,50,50,51:64)]
X.2006.grp<-findInterval(X.2006, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.2006.shad<-shading[X.2006.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2006.shad, exact=T)

#2008 map
X.2008 <- la$pctrep[la$year==2008][c(1:49,50,50,51:64)]
X.2008.grp<-findInterval(X.2008, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.2008.shad<-shading[X.2008.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2008.shad, exact=T)

#2010 map
X.2010 <- la$pctrep[la$year==2010][c(1:49,50,50,51:64)]
X.2010.grp<-findInterval(X.2010, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.2010.shad<-shading[X.2010.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2010.shad, exact=T)

#2014 map
X.2014 <- la$pctrep[la$year==2014][c(1:49,50,50,51:64)]
X.2014.grp<-findInterval(X.2014, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.2014.shad<-shading[X.2014.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2014.shad, exact=T)


par(mfrow=c(3,2),mar=c(0,0,0,0),mai=c(0,0,0,0),oma=c(0,0,0,0),mgp=c(0,0,0),xpd=NA)
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1965.shad, exact=T);text(x=-90,y=32,"1965",cex=2);box()
#map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1970.shad, exact=T);text(x=-90,y=32,"1970",cex=2);box()
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1980.shad, exact=T);text(x=-90,y=32,"1980",cex=2);box()
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.1990.shad, exact=T);text(x=-90,y=32,"1990",cex=2);box()
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2000.shad, exact=T);text(x=-90,y=32,"2000",cex=2);box()
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2010.shad, exact=T);text(x=-90,y=32,"2010",cex=2);box()
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=X.2014.shad, exact=T);text(x=-90,y=32,"2014",cex=2);box()
par(xpd=NA)
legend(x=-96.5,y=32,legend=c("0-5%","5-10%","10-20%","20-30%","30-100%"),fill=shading,cex=1.5)

###OPEN MAP AND ENCODE SPATIAL CONNECTIONS FOR WEIGHT MATRIX###
#call la map
louisiana<-map("county", region=counties.to.map, fill=TRUE, plot=TRUE, resolution=0) 
IDs <- sapply(strsplit(louisiana$names, ","), function(x) x[2])
IDs[50:51]<-"st martin"
la.map.0<-map2SpatialPolygons(louisiana, IDs=IDs)

#create neighbor list from county map
la.nb<-poly2nb(la.map.0)
co.coords<-coordinates(la.map.0) #returns centroids
plot(la.nb, coords=co.coords)

###SPATIAL LAG FOR ALL WAVES USING KRONECKER PRODUCT###
wmat<-nb2mat(la.nb)
N<-64
Time<-49
ident<-diag(1,nrow=Time)
bigW<-(ident%x%wmat)
la$SpatLag<-bigW%*%log(la$pctrep)

###PROBIT VERSION OF TABLE 6.2 WITH OIL###
prob.6.2.2<-lm(qnorm(pctrep)~qnorm(pctrep_t1)+ pctblkreg+pctblack+pctinmig+sla+pcturban+mhhia+lnOil1+as.factor(parish), data=la); summary(prob.6.2.2)
county.clust.2.b<-cluster.vcov(prob.6.2.2,la$parish)
time.fe<-cbind(coeftest(prob.6.2.2,vcov=county.clust.2.b)[1:9,1:3],coeftest(prob.6.2.2,vcov=county.clust.2.b)[1:9,4]/2)
rownames(time.fe)<-c("Intercept","Temporal Lag","Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income ($1000)","Logged Oil Production")
print(xtable(time.fe,digits=4,caption="Parish-Level GOP Registration in Louisiana, 1966-2014. Probit Transform with Fixed Effects and Clustered Standard Errors",align="lrrrr"),type='html')

#FIXED EFFECTS WITH GAS#
gas.6.2.2<-lm(qnorm(pctrep)~qnorm(pctrep_t1)+ pctblkreg+pctblack+pctinmig+sla+pcturban+mhhia+lnGas1+as.factor(parish), data=la); summary(gas.6.2.2)
county.clust.gas<-cluster.vcov(gas.6.2.2,la$parish)
gas.fe<-cbind(coeftest(gas.6.2.2,vcov=county.clust.gas)[1:9,1:3],coeftest(gas.6.2.2,vcov=county.clust.gas)[1:9,4]/2)
rownames(gas.fe)<-c("Intercept","Temporal Lag","Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income ($1000)","Logged Gas Production")
print(xtable(gas.fe,digits=4,caption="Parish-Level GOP Registration in Louisiana, 1966-2014, with Gas Production. Probit Transform with Fixed Effects and Clustered Standard Errors",align="lrrrr"),type='html')


###LOG-LIKELIHOOD OF SPAITAL AND TIME LAG###
spatiotemporal <- function(par, X, y, W, Time, N) {
	n<-N*Time
	rho<-par[1]
	sigma2<-par[2]
	beta<-par[3:length(par)]
	ident<-diag(1,nrow=Time)
	nIdent<-diag(1,nrow=N)
	bigIdent<-ident%x%nIdent
	bigW<-(ident%x%W)
	first<-prod(1-rho*eigen(W)$values)^Time #Accurate and faster to use little W and exponentiate.
    loglikelihood <- log(first)-(n/2)*log(2*pi)-(n/2)*log(sigma2)-((t(y-rho*bigW%*%y-X%*%beta)%*%(y-rho*bigW%*%y-X%*%beta))/(2*sigma2))
    return(loglikelihood)
}

###PROPER SPATIAL MODEL###
la.X<-as.matrix(cbind(1,subset(la,select=c(PROBpctrep_t1, pctblkreg,pctblack,pctinmig,sla,pcturban,mhhia,lnOil1))))

la.spatiotemp <- optim(par=c(0,1,prob.6.2.2$coef[1:9]),	# starting value for beta and sigma
   spatiotemporal,     									# the log-likelihood function
   method="BFGS",               						# optimization method
   hessian=TRUE,                						# return numerical Hessian
   control=list(fnscale=-1,trace=TRUE),    				# maximize function instead of minimize, trace progress
   y=qnorm(la$pctrep), X=la.X,W=wmat, N=64, Time=Time) 	#49 time points
   
#results   
la.spatiotemp$par
sqrt(diag(solve(-la.spatiotemp$hessian)))
la.space.time<-cbind(la.spatiotemp$par,sqrt(diag(solve(-la.spatiotemp$hessian))),la.spatiotemp$par/sqrt(diag(solve(-la.spatiotemp$hessian))),(1-pnorm(abs(la.spatiotemp$par/sqrt(diag(solve(-la.spatiotemp$hessian)))))))
colnames(la.space.time)<-c("Estimate","Std. Err.","z-ratio","p-value")
rownames(la.space.time)<-c("Spatial Lag","Error Variance","Intercept","Temporal Lag","Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income","Logged Oil Production")
print(xtable(la.space.time,caption="Space-Time Model of Louisiana Parish-Level GOP Registration, 1966-2014",align="lrrrr",digits=4),type='html')


###DESCRIPTIVE STATISTICS###
la.X.3<-cbind(la$pctrep,la.X[,-c(1:2)],la$lnGas1)
descriptives<-cbind(apply(la.X.3,2,mean),
	apply(la.X.3,2,sd),
	apply(la.X.3,2,min),
	apply(la.X.3,2,max))
rownames(descriptives)<-c("Proportion GOP Identification","Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income","Logged Oil Production","Logged Gas Production")
colnames(descriptives)<-c("Mean","Std. Dev.","Minimum","Maximum")
print(xtable(descriptives,caption="Descriptive Statistics, Louisiana Parishes, 1966-2014",align="lrrrr",digits=4),type="html")


###UNFOLDING EFFECT OVER TIME: SHOCK###
#1976 to 1977 in St. Tammany Parish
la.X.4<-cbind(1:3136,la$parish,la$year,la.X.3)
la.X.4[la.X.4[,2]=="St. Tammany",1:4]
jump<-qnorm(0.101844862103462)-qnorm(0.0532108284533024)
error.intervention<-rep(0,3136)
error.intervention[756]<-jump

#predictors at mean levels
la.spatiotemp$par
mean.predictors<-matrix(c(1,apply(la.X.3[,-c(1,9)],2,mean,na.rm=T)),ncol=8,nrow=3136,byrow=T)

#matrix algebra
M<-matrix(0,nrow=3136,ncol=3136)
for(i in 1:3072){
	M[64+i,i]<-1
	}
Time<-49
ident<-diag(1,nrow=Time)
bigW<-(ident%x%wmat)
bigI<-(ident%x%diag(1,nrow=64))
rho.hat<-la.spatiotemp$par[1]
phi.hat<-la.spatiotemp$par[4]
beta.hat<-la.spatiotemp$par[c(3,5:11)]
S<-bigI-rho.hat*bigW-phi.hat*M
inv.S<-solve(S)
forecasts.baseline<-pnorm(inv.S%*%(mean.predictors%*%beta.hat))
forecasts.intervention<-pnorm(inv.S%*%(mean.predictors%*%beta.hat+error.intervention))
forecasts.diff<-forecasts.intervention-forecasts.baseline
forecasts.2<-cbind(la$parish,la$year,forecasts.diff)
by(as.numeric(forecasts.2[,3]),INDICES=forecasts.2[,1],FUN=mean)
by(as.numeric(forecasts.2[,3]),INDICES=forecasts.2[,2],FUN=mean)
boxplot(as.numeric(forecasts.2[,3])~forecasts.2[,2])
summary(as.numeric(forecasts.2[,3]))
forecasts.2[forecasts.2[,1]=="St. Tammany",3]
forecasts.2[forecasts.2[,1]=="Washington",3]
forecasts.2[forecasts.2[,1]=="Tangipahoa",3]
forecasts.2[forecasts.2[,1]=="St. Helena",3]
forecasts.2[forecasts.2[,1]=="Livingston",3]


##Map the effect over time##
n.col<-5
br <- c(0,.0001,.001,.005,.01,1)
shading<-gray((n.col-1):0/(n.col-1))

par(mfrow=c(3,2),mar=c(0,0,0,0),mai=c(0,0,0,0),oma=c(0,0,0,0),mgp=c(0,0,0),xpd=NA)
#1977 map
P.1977 <- as.numeric(forecasts.2[forecasts.2[,2]==1977,3])[c(1:49,50,50,51:64)]
P.1977.grp<-findInterval(P.1977, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1977.shad<-shading[P.1977.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1977.shad, exact=T)
text(x=-90,y=32,"1977",cex=2);box()

#1978 map
P.1978 <- as.numeric(forecasts.2[forecasts.2[,2]==1978,3])[c(1:49,50,50,51:64)]
P.1978.grp<-findInterval(P.1978, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1978.shad<-shading[P.1978.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1978.shad, exact=T)
text(x=-90,y=32,"1978",cex=2);box()

#1980 map
P.1980 <- as.numeric(forecasts.2[forecasts.2[,2]==1980,3])[c(1:49,50,50,51:64)]
P.1980.grp<-findInterval(P.1980, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1980.shad<-shading[P.1980.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1980.shad, exact=T)
text(x=-90,y=32,"1980",cex=2);box()

#1990 map
P.1990 <- as.numeric(forecasts.2[forecasts.2[,2]==1990,3])[c(1:49,50,50,51:64)]
P.1990.grp<-findInterval(P.1990, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1990.shad<-shading[P.1990.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1990.shad, exact=T)
text(x=-90,y=32,"1990",cex=2);box()

#2000 map
P.2000 <- as.numeric(forecasts.2[forecasts.2[,2]==2000,3])[c(1:49,50,50,51:64)]
P.2000.grp<-findInterval(P.2000, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.2000.shad<-shading[P.2000.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.2000.shad, exact=T)
text(x=-90,y=32,"2000",cex=2);box()

#2010 map
P.2010 <- as.numeric(forecasts.2[forecasts.2[,2]==2010,3])[c(1:49,50,50,51:64)]
P.2010.grp<-findInterval(P.2010, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.2010.shad<-shading[P.2010.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.2010.shad, exact=T)
text(x=-90,y=32,"2010",cex=2);box()

par(xpd=NA)
legend(x=-96.5,y=32,legend=c("0-0.01%","0.01-0.1%","0.1-0.5%","0.5-1%","1-11%"),fill=shading,cex=1.5)



###UNFOLDING EFFECT OVER TIME: PREDICTOR###
#1968 to 1969 in West Feliciana Parish
la.X.4[la.X.4[,2]=="West Feliciana",]
jump.2<-0.573356568813324-0.287450969219208
intervention.predictors<-mean.predictors<-matrix(c(1,apply(la.X.3[,-c(1,9)],2,mean,na.rm=T)),ncol=8,nrow=3136,byrow=T)
intervention.predictors[255,2]<-intervention.predictors[255,2]+jump.2

#matrix algebra
M<-matrix(0,nrow=3136,ncol=3136)
for(i in 1:3072){
	M[64+i,i]<-1
	}
Time<-49
ident<-diag(1,nrow=Time)
bigW<-(ident%x%wmat)
bigI<-(ident%x%diag(1,nrow=64))
rho.hat<-la.spatiotemp$par[1]
phi.hat<-la.spatiotemp$par[4]
beta.hat<-la.spatiotemp$par[c(3,5:11)]
S<-bigI-rho.hat*bigW-phi.hat*M
inv.S<-solve(S)
forecasts.baseline<-pnorm(inv.S%*%(mean.predictors%*%beta.hat))
forecasts.intervention<-pnorm(inv.S%*%(intervention.predictors%*%beta.hat))
forecasts.diff<-forecasts.intervention-forecasts.baseline
forecasts.X.2<-cbind(la$parish,la$year,forecasts.diff)
by(as.numeric(forecasts.X.2[,3]),INDICES=forecasts.X.2[,1],FUN=mean)
by(as.numeric(forecasts.X.2[,3]),INDICES=forecasts.X.2[,2],FUN=mean)
boxplot(as.numeric(forecasts.X.2[,3])~forecasts.X.2[,2])
forecasts.X.2[forecasts.X.2[,1]=="West Feliciana",3]
forecasts.X.2[forecasts.X.2[,1]=="East Feliciana",3]
forecasts.X.2[forecasts.X.2[,1]=="Concordia",3]
forecasts.X.2[forecasts.X.2[,1]=="Pointe Coupee",3]
forecasts.X.2[forecasts.X.2[,1]=="East Baton Rouge",3]
forecasts.X.2[forecasts.X.2[,1]=="West Baton Rouge",3]
summary(as.numeric(forecasts.X.2[,3]))

##Map the effect over time##
n.col<-5
br <- c(0,.0001,.001,.0025,.005,.01)
shading<-gray((n.col-1):0/(n.col-1))

par(mfrow=c(3,2),mar=c(0,0,0,0),mai=c(0,0,0,0),oma=c(0,0,0,0),mgp=c(0,0,0),xpd=NA)

#1969 map
P.1969 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==1969,3])[c(1:49,50,50,51:64)]
P.1969.grp<-findInterval(P.1969, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1969.shad<-shading[P.1969.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1969.shad, exact=T)
text(x=-90,y=32,"1969",cex=2);box()

#1973 map
P.1973 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==1973,3])[c(1:49,50,50,51:64)]
P.1973.grp<-findInterval(P.1973, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1973.shad<-shading[P.1973.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1973.shad, exact=T)
text(x=-90,y=32,"1973",cex=2);box()

#1980 map
P.1980 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==1980,3])[c(1:49,50,50,51:64)]
P.1980.grp<-findInterval(P.1980, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1980.shad<-shading[P.1980.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1980.shad, exact=T)
text(x=-90,y=32,"1980",cex=2);box()

#1990 map
P.1990 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==1990,3])[c(1:49,50,50,51:64)]
P.1990.grp<-findInterval(P.1990, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.1990.shad<-shading[P.1990.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.1990.shad, exact=T)
text(x=-90,y=32,"1990",cex=2);box()

#2000 map
P.2000 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==2000,3])[c(1:49,50,50,51:64)]
P.2000.grp<-findInterval(P.2000, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.2000.shad<-shading[P.2000.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.2000.shad, exact=T)
text(x=-90,y=32,"2000",cex=2);box()

#2010 map
P.2010 <- as.numeric(forecasts.X.2[forecasts.X.2[,2]==2010,3])[c(1:49,50,50,51:64)]
P.2010.grp<-findInterval(P.2010, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
P.2010.shad<-shading[P.2010.grp]
map(database="county", region=counties.to.map, fill=T, plot=T, interior=T, col=P.2010.shad, exact=T)
text(x=-90,y=32,"2010",cex=2);box()

par(xpd=NA)
legend(x=-96.5,y=32,legend=c("0-0.01%","0.01-0.1%","0.1-0.25%","0.25-0.5%","0.5-1%"),fill=shading,cex=1.5)




###SPATIAL MODEL WITH GAS###
la.gas<-as.matrix(cbind(1,subset(la,select=c(PROBpctrep_t1, pctblkreg,pctblack,pctinmig,sla,pcturban,mhhia,lnGas1))))

gas.spatiotemp <- optim(par=c(0,1,prob.6.2.2$coef[1:9]),	# starting value for beta and sigma
   spatiotemporal,     									# the log-likelihood function
   method="BFGS",               						# optimization method
   hessian=TRUE,                						# return numerical Hessian
   control=list(fnscale=-1,trace=TRUE),    				# maximize function instead of minimize, trace progress
   y=qnorm(la$pctrep), X=la.X,W=wmat, N=64, Time=Time) 	#49 time points
   
#results   
gas.spatiotemp$par
sqrt(diag(solve(-gas.spatiotemp$hessian)))
gas.space.time<-cbind(gas.spatiotemp$par,sqrt(diag(solve(-gas.spatiotemp$hessian))),gas.spatiotemp$par/sqrt(diag(solve(-gas.spatiotemp$hessian))),(1-pnorm(abs(gas.spatiotemp$par/sqrt(diag(solve(-gas.spatiotemp$hessian)))))))
colnames(gas.space.time)<-c("Estimate","Std. Err.","z-ratio","p-value")
rownames(gas.space.time)<-c("Spatial Lag","Error Variance","Intercept","Temporal Lag","Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income","Logged Gas Production")
print(xtable(gas.space.time,caption="Space-Time Model of Louisiana Parish-Level GOP Registration, 1966-2014, with Gas Production",align="lrrrr",digits=4),type='html')

